{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Formate perovskites\n",
    "\n",
    "Perovskites are an important class of inorganic materials, particularly for solar energy applications.\n",
    "Usually, the perovskite structure contains an A cation, a B cation and an C anion in the ratio 1:1:3 (black, blue and red in the figure below, respectively, [see wikipeidia for more information](https://en.wikipedia.org/wiki/Perovskite_(structure)).\n",
    "\n",
    "<img src=\"Perovskite.jpg\">\n",
    "\n",
    "\n",
    "Here we search for charge inverted perovskites, i.e. with an anion on the A site. This class of material is closely related to perovskites, and may represent another fruitful search space for new photovoltaic materials. \n",
    "\n",
    "In this example  we assume a simple [formate moelcule](https://en.wikipedia.org/wiki/Formate) as the C-site and uses [Goldschmidt ratio rules](https://en.wikipedia.org/wiki/Goldschmidt_tolerance_factor) as part of the screening.\n",
    "These rules allow us to estimate whether or not a perovskite structure is likely to form based on data about ionic size alone. The tolerance factor is defined as a ratio of the radii of the A, B and C species\n",
    "\n",
    "\n",
    "\\begin{equation*}\n",
    "t = \\frac{r^A + r^C}{\\sqrt{2}(r^B + r^C)} ,\n",
    "\\end{equation*}\n",
    "\n",
    "Values of t > 1 imply a relatively large A site favoring a hexagonal structure, 0.9 < t < 1 predicts a cubic structure, and 0.7 < t < 0.9 means that the A site is small, preferring an orthorhombic structure. For t < 0.7, other (non-perovskite) structures are predicted.\n",
    "\n",
    "We also apply the standard charge neutrality and electronegativity tests as described [in the docs](https://smact.readthedocs.io/en/latest/examples.html#neutral-combinations).\n",
    "\n",
    "\n",
    "We outline 2 approaches for achieveing the same result:\n",
    "\n",
    "1. The methodology is written out explicitly for transparency, including accessing the `smact` data directory directly and storing element and species information as simple lists of strings, ints and floats.\n",
    "\n",
    "2. Fewer lines of code are used, making use of inbuilt `smact` functions. Element and species information is stored directly as Element and Species objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "### IMPORTS ###\n",
    "get_ipython().magic(u'matplotlib inline') # plots appear in the notebook\n",
    "import smact\n",
    "from smact import Species, Element\n",
    "import smact.lattice as lattice\n",
    "import smact.screening as screening\n",
    "\n",
    "from itertools import product\n",
    "import copy\n",
    "import os\n",
    "import csv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from smact import data_directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we define the positions of the 3 sites in the perovskite structure and specify the allowed oxidation states at each site. Note that the A site is defined as an anion (i.e. with a -1 oxidation state)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "site_A = lattice.Site([0,0,0],[-1])\n",
    "site_B = lattice.Site([0.5,0.5,0.5],[+5,+4])\n",
    "site_C = lattice.Site([0.5,0.5,0.5],[-2,-1])\n",
    "perovskite = lattice.Lattice([site_A,site_B,site_C],space_group=221)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Approach 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now search through the elements of interest (Li-Fr) and find those that are allowed on each site. In this example, we use the F- anion with an increased Shannon radius to simulate the formate anion. We access the Shannon radii data directly from the smact data directory and are interested in the octahedral (6_n) Shannon radius. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "search = smact.ordered_elements(3,87) # Li - Fr\n",
    "\n",
    "A_list = []               # will be populated with anions\n",
    "B_list = []               # will be populated with cations\n",
    "C_list = [['F',-1,4.47]]  # is always the \"formate anion\"\n",
    "for element in search:\n",
    "    with open(os.path.join(data_directory, 'shannon_radii.csv'),'r') as f:\n",
    "        reader = csv.reader(f)\n",
    "        r_shannon=False\n",
    "        for row in reader:\n",
    "            if row[2] ==\"6_n\" and row[0]==element and int(row[1]) in site_A.oxidation_states:\n",
    "                A_list.append([row[0],int(row[1]),float(row[4])])\n",
    "            if row[2]==\"6_n\" and row[0]==element and int(row[1]) in site_B.oxidation_states:\n",
    "                B_list.append([row[0],int(row[1]),float(row[4])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*NB: We access the data directly from the data directory file here for transparency. However, reading the file multiple times would slow down the code if we were looping over many (perhaps millions to billions) of compositions. As such, reading all the data in once into a dictionary, then accessing that dictionary from within a loop, could be preferable, e.g.:*\n",
    "``` python\n",
    "for element in search:\n",
    "    ...\n",
    "    r_shannon = shannon_radii[element][coordination]\n",
    "    ...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We go through and apply the electronegativity order test (pauling_test) to each combo. Then, we use Goldschmidt tolernace factor to group into crystal structure types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We define the different categories of list we will populate\n",
    "charge_balanced = []\n",
    "goldschmidt_cubic = []\n",
    "goldschmidt_ortho = []\n",
    "a_too_large = []\n",
    "A_B_similar = []\n",
    "pauling_perov = []\n",
    "anion_stats = []\n",
    "\n",
    "# We recursively search all ABC combinations using nested for loops\n",
    "for C in C_list:\n",
    "    anion_hex = 0\n",
    "    anion_cub = 0\n",
    "    anion_ort = 0\n",
    "    for B in B_list:\n",
    "        for A in A_list:\n",
    "            # We check that we have 3 different elements\n",
    "            if B[0] != A[0]:        \n",
    "                # Check for charge neutrality\n",
    "                if int(A[1])+int(B[1])+3*int(C[1]) == 0:\n",
    "                    charge_balanced.append([A[0],B[0],C[0]])\n",
    "                    # We apply the pauling electronegativity test\n",
    "                    paul_a = smact.Element(A[0]).pauling_eneg\n",
    "                    paul_b = smact.Element(B[0]).pauling_eneg\n",
    "                    paul_c = smact.Element(C[0]).pauling_eneg\n",
    "                    electroneg_makes_sense = screening.pauling_test([A[1],B[1],C[1]], [paul_a,paul_b,paul_c])\n",
    "                    if electroneg_makes_sense:\n",
    "                        pauling_perov.append([A[0],B[0],C[0]])\n",
    "                        # We calculate the Goldschmidt tolerance factor\n",
    "                        tol = (float(A[2]) + C[2])/(np.sqrt(2)*(float(B[2])+C[2]))\n",
    "                        if tol > 1.0:\n",
    "                            a_too_large.append([A[0],B[0],C[0]])\n",
    "                            anion_hex = anion_hex+1\n",
    "                        if tol > 0.9 and tol <= 1.0:\n",
    "                            goldschmidt_cubic.append([A[0],B[0],C[0]])\n",
    "                            anion_cub = anion_cub + 1\n",
    "                        if tol >= 0.71 and tol < 0.9:\n",
    "                            goldschmidt_ortho.append([A[0],B[0],C[0]])\n",
    "                            anion_ort = anion_ort + 1\n",
    "                        if tol < 0.71:\n",
    "                            A_B_similar.append([A[0],B[0],C[0]])\n",
    "anion_stats.append([anion_hex,anion_cub,anion_ort])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 40, 91]]\n",
      "Number of possible charge neutral perovskites from Li to Fr = 132\n",
      "Number of Pauling senseibe perovskites from Li to Fr = 132\n",
      "Number of possible cubic perovskites from Li to Fr = 40\n",
      "Number of possible ortho perovskites from Li to Fr = 91\n",
      "Number of possible hexagonal perovskites from Li to Fr = 1\n",
      "Number of possible non-perovskites from Li to Fr = 0\n",
      "----------------------------------------------------------------\n",
      "Structures identified with cubic tolerance factor 0.9 < t < 1.0 \n",
      "----------------------------------------------------------------\n",
      "Cl C (HCOO)3\n",
      "Br C (HCOO)3\n",
      "Cl Si (HCOO)3\n",
      "Br Si (HCOO)3\n",
      "I Si (HCOO)3\n",
      "Cl S (HCOO)3\n",
      "Br S (HCOO)3\n",
      "I S (HCOO)3\n",
      "I Ti (HCOO)3\n",
      "Br V (HCOO)3\n",
      "I V (HCOO)3\n",
      "Br Cr (HCOO)3\n",
      "I Cr (HCOO)3\n",
      "Br Mn (HCOO)3\n",
      "I Mn (HCOO)3\n",
      "I Fe (HCOO)3\n",
      "Br Co (HCOO)3\n",
      "I Co (HCOO)3\n",
      "Br Ni (HCOO)3\n",
      "I Ni (HCOO)3\n",
      "Br Ge (HCOO)3\n",
      "I Ge (HCOO)3\n",
      "Br Se (HCOO)3\n",
      "I Se (HCOO)3\n",
      "I Zr (HCOO)3\n",
      "I Nb (HCOO)3\n",
      "I Mo (HCOO)3\n",
      "I Tc (HCOO)3\n",
      "I Ru (HCOO)3\n",
      "I Rh (HCOO)3\n",
      "I Pd (HCOO)3\n",
      "I Sn (HCOO)3\n",
      "I Tb (HCOO)3\n",
      "I Hf (HCOO)3\n",
      "I Ta (HCOO)3\n",
      "I W (HCOO)3\n",
      "I Re (HCOO)3\n",
      "I Os (HCOO)3\n",
      "I Ir (HCOO)3\n",
      "I Pt (HCOO)3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAD1CAYAAABnVo9yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8W9d99/HPwSAFiqJEiVrWljWtYW1ZtqanZDmyHduxk9YjqVO6SZMUbp88KdqmSZOgTtuHt0ma4TQeGXY8Ysfbkpcka+9ha1OyNklRi1ocAO55/riXexPjAuDv/XrxBRG4uPiRhL48PPcMpbVGCCGE81xOFyCEEMIigSyEEElCAlkIIZKEBLIQQiQJCWQhhEgSEshCCJEkJJBFs5RSh5VSWik1v5XjVtjHPZKYyoRITxLIQgiRJCSQhRAiSUggCyFEkpBAFnGllJqplHpBKXVcKVWllCpVSr2hlJrdxLFP2n3RS5VSqonHn7Uff7Opx4VIdRLIIm6UUn8PrAO+ABQDrwOFwGJgpVLqqw2e8i1gB3Ab8J0G53oYeBg4BjysZREWkYYkkEVcKKUWAv8FFAGztNbTtNb3aa1nAXOBS8DPlVKjqp+jta7ACu9LwL8ppW6wzzUW+AUQBr6otT6b2K9GiMSQQBZtsdzuKmjyA5jXxHO+b98+qrXeUPcBrfUa4AeAF8hv8Nh++z4P8IJSahDwMpAF/LP9XCHSksfpAkRKWIbV5dCchUDf6k+UUnnAdOAC8F4zz1lp385q+IDW+nl77PNXgU+A7sBS4D/aW7gQqUQCWbTFE1rrFc09qJRaQZ1ABoYBCsgBwq1cf+vdzP3fBG4FhgCngAel31ikOwlkEQ9u+7YMeK2VY083c/9cYLD9757AyBaOFSItSCCLeDhm34a01o+098lKqf7A77Fa2c8AX8bqT56ktT4XsyqFSDJyUU/EnNb6BFbfb15r62A0pJRyAc8BfYCfaK2/AvwOq7X8TIxLFSKpSCCLePkX+/YPSqlbGz6olMpQSi1RSjW8qPddYAGwGfi2fd/XgL3AnUqpb8WrYCGcJoEs4kJr/Trw90A/YJlSap89Q+9PSqkNWBfqXgeurX6OUmoBVpBfAB7QWlfZ57qMNT65AvgPpdS0xH41QiSGBLKIG611ATAVeArrQt8tWLPwcrGGvX0VeAlAKdUHq6vCBXxVa32wwbk+Af4OyABeVEp1T9CXIUTCKBlJJIQQyUFayEIIkSQkkIUQIknIOGSRMgryfW6gh/dAVk5mYdccwIfVp+w6d/PNZqR7d4AIEMK6MHgWOBcwjJBTNQvRHtKHLJJGQb4vG7gGGAeMB8ZijdLIxZqt1w1QrrPej7M29Jhb97lnlizZpb3ecc2c+hJwDiugS4CDDT4OBQzjcuy/IiHaR1rIIuEK8n1Z1AZv3Y/BWLPzWuYx29vVlm1/DGrugKDfX4y1FvOm6o+AYRS183WEiIoEsoi7gnxfN6whb4uB+dQuPtQh2q3j8b7tZ3/cVn1H0O8/iRXOa7FWrdsRMAz5k1LEjXRZiLgoyPeNxArgxVgLBWXE6tyqwrW56/Je9SaHtNJlESslWMG8DHgvYBilcX490clIIIuYKMj3ebGCdzFwB9bqbPFRpXZmf5g3se5dCQrkujSwBWvx/D8GDONYK8cL0SoJZNFh9qiHzwEPYnVJdEvIC4fZl/1+79F173IgkOvSwGrgeeDlgGGccagOkeIkkEW7FeT7crGmPX8NawH5xDI5kr2sd73XdTiQ6wphdWn8L/BWwDBMh+sRKUQCWbRZQb5vAvAN4C+w9rhzhqY0e2nvejuNJFEg13UI+DnwVMAwypwuRiQ/CWTRIrtb4k6sIJ7vbDU2zeXspb271r0rSQO52mWsNZ1/FjCMPU4XI5KXBLJoUkG+rye13RKDWzk8sTS669I8FLWb9SV5IFfTWEuOfjdgGJ84XYxIPhLIoh57zPC3gcdxsluiFV3f63VZRVw1reQUCeRqGmvZ0e8FDGOv08WI5CGBLAAoyPd5sFrE38PaPimpZX3U85Sr0l1TZ4oFcrUI1siM7wcM42BrB4v0J6u9CQryfXcCnwK/IAXCGACvLne6hBhwYw0Z3Bv0+38S9Pt7OF2QcJa0kDuxgnzfKOBnQKM975Kdb12Pfe7z3pqxyCnaQm6oFPgn4DcyRbtzkhZyJ1SQ7/MV5Pt+iLUzdMqFMYD2mpVO1xAHvYFfA2uDfv+1rR0s0o8EcidTkO9bAuzBaonFbH2JRNNeawPUNHUdsCXo9/9X0O/v4nQxInEkkDuJgnxfVkG+72msYVeJn10XY9qr033ReTfWrt1bgn7/ZKeLEYkhgdwJFOT7xgIbgS87XUusaK8ZcbqGBLkG2BD0+/8x6PfL/9c0Jz/gNFeQ73sIa03fVL/gVZ9Xd6Y1IrxAEFgZ9PuHOV2MiB8J5DRld1E8A/wW6Nra8alGe83OOAphNrAj6Pd/welCRHxIIKehOl0UjzhcStxoT2fMY8Ba4vTFoN//Y+nCSD/yA00zBfm+h0nHLooGtLfTBnK1bwNvB/3+XKcLEbEjgZwm7C6KZ4FnScMuikY8ptvpEpLAQmBj0O9P61++nYkEchooyPd1B94HHna6lkSJ00anqWgEsD7o93/O6UJE9CSQU1xBvi8PWA5c73QtCeWWHdPryAZeDfr9neYXcrqSQE5hBfm+/sBKoNNNHNAunbKzDOPEAzwT9Psfd7oQ0XESyCmqIN83BFiFNXGg83HpTKdLSEIK+H9Bv//fnS5EdIwEcgoqyPeNxArjq52uxTEufE6XkMS+E/T7fy3D4lKP/MBSTEG+bzzwMTDI6VocpSSQW/FV4Kmg369aPVIkDbkwkkIK8n1TsbaY7+V0Le2xtTDC2t1his5oTA19eiimj3Zz/TVuXKp9eXGlUrNiR5hPj5i9zp0/jqk1OW43/d59d+CC667jqj6N19e/ePkyb69cyYHDhzG1ZtjAgSyeN49ePZpeD/6tFSvYuHMnX//Sl+ibl9ehrzlJPAJcBL7pcB2ijRxtISulbldK/V4pVaiUuqSUqlBKHVdKva2Uekwp1S0Gr/GsUkorpR5p5/O+Zz/ve9HWEAsF+b4bgI9IsTB+ZXWI5z4KcaxUM6y/i1EDXZSWaf68Jszv3g9htmMG9LlLmoJXKvlwe4RL5do1KjOTiT4fbqXYfehQ91++8AKfHjhQ7zlaa373+ut8sn8/fXr1YlD//uz77DN+86c/UV5R0eg1TpSUsGHHDmZPnZrqYVztG0G//0dOFyHaxpEWslKqD9Ymj/Psu/YA7wFVwEDgZuB24AdKqWla6yNO1Jks7KnQ7wA5TtfSHjsPRVi7O0I3H3x9SQa9u1u//y9e0fzyrSo+OWyyeleEuRPa9jZ8e0OIc5dg7CAXD93iJfejXpdUxJVtas1zeXmlq7dv7/36hx8ydvhw3G5r3sjugwc5eeoU08aP5+6bbwbgw/Xr+Wj9ejZ/+ilzpk2rOb9pmrz24Yf0yMlhwcyZMf5uOCoQ9PsvBgzjCacLES1LeAtZKdUDWIMVxuuAa7XW12itP6+1fkBrPRvIA/4v1gLqTk0N/R9grH3rmIJ8X0/gTVIsjAE+3B4G4I6Z3powBuiWpbhntheAj7aHMdu4jVjhSWuBt1umeMjwKLRHXwFwKcXCWbNOeT0erlRUcOb8+ZrnnDx1CoCp19QORpk+fjwAR4uK6p1/3fbtnDx1ijtvugmvJ+168/496Pd/zekiRMuc6LL4H6zZRRuBG7XWOxseoLW+qLX+D2AqUJLg+qprOK213qu1Pu3E60PNTtAvk4KjKc5f0hw/rXG74Nrhjd9mV1/lontXuFgOR0raFsiehpOlm9noNMtXe73vit0t4etSu/FG9b/Dkdollc9fvMgH69YxaexYRgwe3KZ6UtDPgn7/7U4XIZqX0EBWSl0NfNH+9DGtdeNOvDq01oVa6yL7uSvsPt35zZy71b5ipdQkpdRrSqnTSqkrSqktSqkmF21vrQ9ZKTVWKfVru/+7XCl1Tim1Uyn1X0qpWO3I8RPgxhidK6FOnLFas/1yFV5P0xfuBvV21Tu2NWMGWYn8/rYwVWGN9ljvH601S9et6xMKhxkzfDjZWVk1z8nNsf6wKD17tua+6n9XPwbw5vLleNxubp8zp21fYGpyAX8M+v2dc+x6Ckh0C/kO+zU/0VpvS/Brz8TqIhmPte7DOuBa4Gml1E/bcyKl1EPAdqyhRQp4C2vGnAtr250F0RZbkO97DEjZPzHPXrRavbndmh9FkZut6h3bmoXTPQzurdhz1OSHz1fy35+d6f/k6dN8t7iYtTt35k0aM4YvLFxY7zmjhw1DKcWH69dz5vx5Lly+zNJVqwAYO3w4ALsKC9l76BCL5syha50wT1M5wBtBvz+lLg53FonuKJtq325K8OsCPAb8FHhcax0BUErNxArnbyillmqt32ntJEqp6cBTWEH8KPC01rWdoEqpsdEWWpDvWwD8LNrzOKnS3vEuo4V3WPVjlW3crjS7i+KxOzJ4dU2IzftN9hPKAeuFcnNyQkMHDszIzKg/o7pvr17cMHkyq7dupeDZZ2vunzBqFCOHDqWyqoq3Vqxg2MCBTBlXu2iaaZqYpokn/fqSweoC+1PQ7781YBjpvjdhSkn0u623fXsqwa8LcBL4dnUYA2itNyilDOC7gB9rJENr/gnr+/ZjrfVTDR/UWu+JpsiCfN9wrH7j1E4C+1dULGcllJw3eXpZiMoqzZcWeBl9sevm7kVdpx2pquK5SMR87YMPOHryJPfcemu95y2aO5fhgwZx4MgRTNNk2MCBjB85EoD31qzhSnk5f3XPPQCcOX+eN5cv5+CxY5imSd9evVg0Zw4jhw6N4VeSFOZjNVD+xuE6RB2daabey1rryibu/719O1sp1WIIKqXcWEPyAH4Ty+IACvJ93bBGVKT8n5OZdkO1Mtz8MVXh+se2JGJqfvt+iDNlmkduzWDqSDe53VwVPpeLMV268Ohddx3Ozspi6+7dHDp2rNHzRw8bxh3z57PkxhuZMGoUSimOFxezYedO5s2YQV5uLhWVlTz1yiscLy7mjvnz+eLixSil+N0bb3C8xJFry/H2WNDv/wunixC1Eh3IpfZt4+lU8fdZM/cfBUygC60HYR7W4u9hrXVhDGujIN/nAv5ImiwWVN0/fK6F/uHzl6zHema33o4+ekpTck7Ts5tiaF/7bVtno9Nsny8yym7FFh492ur5qscc5+XmMtcei7xj3z7KLl7ktjlzmDlxIuNHjuSLixdjmiZrtmxp9Zwp6heycWrySHQgV7+rp8fh3LH4WpzcF+jvgcUOvn5MDcizfhzF5zShcNPf1qOlZr1jW3LODu8udVrTDTc67ZJpLQDX1Ay8htZu20ZxaSl33XQTHnsSSXGp1V4Y3L9/zXF5ublk+XwUnXZs9GO85QDPB/3+1O4iSxOJDuS3sVqjE5RS7V3Dt/rST3Yzj7c21GxoM/cPxvo+VABnmzmm2mngCuCxh/DFhD0T799idb5kkJutGJiniJiw41DjYW0HT5qUXYZuPhjSt/UWcvcs65hT5zXllVYOa2/9GSXH7Ikeud27t3iu8xcu8OH69UwdP56hAwbU3J/htSarhEK117m01oRCoZj2hSeh64B/dboIkeBAtv/Mf9H+9JdKqRbXtFVKXa2Uqm6unLBvxzRxXF9gSisvf59Sqqneyuo+tDVa6xZ6PMG+IPiB/emjrbxemxTk+9xY++B1aeXQlHPjJKvR9daGEKfLakP5YrnmlTWhmmPqLjC0+tMwT7xYyfPL6w+9GNJXkZMFoQi8uDJERZVGe7QLwLTGIfc+VlyMy+Vi3IgRLdb1xvLlZHi9LJw9u979/ey1K7bs3l1z3yf79xMKh+nfxKJFaSYQ9PvnOl1EZ+fEnyl/izUmeCbwkVLqMa31J3UPUEp1BfKxRj/MB4qAD7FWr/q6Uuq5OhNGegK/pfmWc7UBwBNKqX/Q2up7tIewVe+w8JM21v8jrK6Ff1BK7dNaP9ug9jEAWuu9bTzft4EZbTw2pVw73M3115is3R3hP/9UxagBLtwuOHDCpCIE44e6mD2u/vS7yxWa0jJNToPhwB634oH5Xp5ZFuKTwyYH/1hJ326hCd0vhzleVcXp48f7KKVaXMUN4NMDB9j32Wfcv2hRvdl7ABNGj2bFpk1s3LmT4tJSfF26UHjkCB63m9lTWvt9n/JcwO+Dfv+4gGFccrqYzirhgay1PquUmo21uNBsYKdSajewF6tbYgBWQGViTZuu7kZ4CSs8JwO7lFJrsNa6mI41pO014K4WXvpXWBMtPqeU2ow1BG8e1vfgF1rrN9tY/0al1F8DTwLPKKX+GatvPBNrfOd44Mv219OignzfOOB7bXndVHXPbC/D+rlYsyvMwSITHcXym6MHuvn7exUrd0YoPGly9KzZHbOcHLebsUOHls2dObN73f7fhioqK3lrxQpGDhnCxNGjGz3ucbt55O67WbpqFYVHjhCORBjYrx+33nBDk8t6pqHBwPexrmcIByjdxoVd4vLiSt2BNZV6FtAXKxxLsWbBvQ48r7W+XOf4nsAPgTuxArUIeBWr/+unWLsuf7luq1Up9Wz1/cBOrDfcDYAPa5W5XwBP6QbfCHvK9L8C39daf6+J2idg/YK4EegHXMYasfEe8D9a6xYv9Rfk+xSwms62OWkMqQrX5q7Le00DOLNkyS7t9Y5r7TmiVRFgWsAwtjtdSGfkaCB3ZgX5vkeB/3W6jpQWUjuzP8ibCBLIMbYRmBUwjLYtMiJipjNNDEkaBfm+PODHTteR8pROuwuhSWIGMoPPERLIzvhPoKfTRaQ8V/qNTEkiwaDf33yHvIgLCeQEK8j3zcbq0xbRko1O4ykHCDpdRGcjgZx4/05s19zpzNJ+rUyHPRT0+yc4XURnIoGcQAX5vvlYQ/1EbGRp5Kp0HLmQax0JJYGcWP/sdAFpRaFw1w6LFHGxKOj3p/U2KslEAjlBCvJ91wE3OV1Huqne6FTE1Q+dLqCzkEBOnH9xuoC01MxGpyKm5gb9/lucLqIzkEBOgIJ832RAdvuNg+qNTkXc/aPTBXQGsgZqYjTbd7zq0zCfFZsUndVcKtdUVIEvE67q5WL6KDdTRrhQTaz3YGrN2t0RNu2LcOq8xqWgfy/F9dd4mDLC3cQrtayjdRw9ZfLG+hDHSjVdvDBxuJs7ZnrI9DZRs6n5yWtVlFfC/7kvo9ndqNtDe00J5MRYEPT7J8mU6viSQI4zewGhu5t7fPmOMJfKoV+utRNGhsdajL3whMmBEyY7Drl45FZvvUV4TFPz7Pshdh0x6eKF0QNdhCNw4KTJcx+FOFJicvcN3nbV2ZE6yi5rfvlWFRHTquHcJeuXxNmLmq8uarzS6apPIxw/rfnr270xCWMA7dWySWfiPA485HQR6UwCOf7+iRbGHf/lTRkM6KUatSiLz5r86u0qdh0x2bw/wozRtT+qjz+NsOuISd9cxd8szqCbvXh7aZnJz9+oYvWuCCMHuBg/tO0t5Y7UsXxHmKowfGmBl6kj3Zim5tfvhth7zOToKZPBfWp7xM5d0izbEmbKCBejB7a/Bd8c7TUlkBPngaDf/52AYZx0upB0JX3IcVSQ7xsJfKGlY4b3czX5532/ni5uGGeF3/7jtWu8mKZm+Q5rHf17Zntrwhigd3cXi2daLeMPtrW41n7UdQAcP23iccPkEdbbyOVSzBhthe2RU/WP/fOaEG4XLJnVvpZ7a7S3dhdxEXde4BtOF5HOJJDjyw90uDnosvPR464NysOnNJfKoXtXuLp/4x/ftcOtReCPlWrKLsdmzkRTdQBcqQBfBvW6MbLsPWBCdX4f7PzMatHfMdNDN1+MJynW2ehUJER+0O/v6nQR6UoCOU4K8n0eWmkdt+TMBZN1e6xUGzek9sd04rSVP4N7N/2jy/Ao+uZaoXfiTPRZ1VwdALndFJcq4EpFbfCfOm/9u1eOVUNFlea1tSGG969tPcdSw41ORdzlAvc5XUS6kj7k+LkF6NXWgzfus3bUiJjWxbLDxRoN3DTJzYRhtUF29qKVP7nZzbc0c7MVJ89ozl5of1a1tQ6wAnrvMZPX1oa483ov5y5pVuwMk+mtbb2/u8m6WJi/2NvkKI1oNdzoVCTEQ1j7QIoYk0COny+25+DPijWb99e2aF0KFk7zMG9C/RCssi9hZbTQFVv9WGUHLne1tQ6AmWPcbN4fYUuhyZbCypr775ntIdunOFZqsmZ3hFumeOjbo7Z1HY5Yw/RcrhgMe7M3OhUJNT/o9w8OGEaLu+KI9pNAjoOCfF8XWt7fr5H753m5f56XUFhz5qJm074I720Js+NQhEcXZtC9qxVe8W4OtrUOALdL8bXPZbD5QITjpZoML0wc5mZoXxemqXn54xB5OYqbJllhvvdYhLc2hCk6awXyiAEu7rnBQ173KDLVY0ogJ54CHsTa8FfEkLyZ42Mx0K0jT/R6FP1yXXzuOi+3z/Bw8ozmz2tqm7qZduu3qoXWb/VjmVEMaGitjmoet+K6MR7uneNlyXVehva13lIffxrh5BnNfXO8eNyKo6dMnlpqPf+hm718fraH46Umv3yrioqqjv+a0W4tjQpnPOh0AelIAjk+2tVd0ZwZo6yW5a4jJhH72lXPblYL9dyl5kPs/KX6x8ajjpacu6RZtjnM9NFurr7Keout3BnG1PDILV6uHe5m1lgPi2d6OH8ZthVGMXLNLX/lOWR00O+f6XQR6UYCOcYK8n3dsFrIUeuSafXhmtoaYgYwIM/6kR0tbXoERVVYU3xO1zs2HnW05NXVITK88LmZtVl58qymaxfqdU8MtSeOnDwbVQu58ZRAkSj3O11AupFAjr27IDZ7vR0qMjG1Nda3q33GoX0U2T4ouwwHixqH8o5D1giJQb1Vvf7eWNfRnB2HIuw+anLnLC9ZXWpfP8NTf2wyQFX75q40TTY6dVJMGh6ilgRy7D3Q1gMPFZlsORAhHGncQvys2OSllVZizRjjrhmR4HIp5k+0Wp6vrA5xsbz2uaVlJm9vsPppb57c+C/5tzeGeOLFSt7eWL8vuCN1NKV6zPGoga5GCxz17+miKgzbD9Z2T2zab/17QF4Uvzhko1MnjQr6/SOcLiKdSP9bDBXk+3phjT9uk9MXTF5cGebPa6zuhW5ZUFkFZy5qSuxuh7GDXSyaVv/HNG+Cm0NFJruPmjzxQiUjB7iImLD/hEk4ArPHuZtcx+LCFU1pmebClfrB29E6GnpnY5jySrh3duPjbprkZmthhOc+CrH5QITySs3hEk2vborJV0cxYUQ2OnXaYuAnTheRLiSQY+sOrPn+bXJ1fxe3THFzqEhTWmZyuMS6v5sPJg5zMWWkmwlNBKvLpfjybV7W7oqwaX+EfcdNlIKBeYobxrV/+c2O1lHX0VMma/dEWDTdQ6+cxn949e7hIv/2DN7ZFOLACROPy5rm/bnrvE2uodEOstGpsySQY0hpmegUMwX5vl8DX3W6jk5Fo7suzePskjt3a693nNPldEJVQK+AYVxyupB0IH3IsSXDgBJNNjp1WgZwo9NFpAsJ5BgpyPdlA9JCc4BsdOq4G5wuIF1IIMfOdKJYalNEQTY6ddp1TheQLiSQY0felA6RjU4dNy3o98sAgRiQQI4d6T92iGx06rgsYKLTRaQDCeTYkUB2iGx0mhRmOV1AOpBAjoGCfN9QoJ/TdXRWstFpUpAuuxiQQI4NeTM6SDY6TQqTnC4gHUggx4Z0VzhJNjpNBiOCfr/kSZTkGxgbEsgO0l5TAtl5XYDBTheR6iSQY2Ok0wV0ZtqrY797quiIUU4XkOokkKNUkO/zAXlO19GZaY8EcpIY7XQBqU4COXqDnC6g05ONTpOFtJCjJG/k6Em/mcNko9OkIV13UZJAjp60kJ0mG50mCxmLHyUJ5OgNdLqAzk67dabTNQgAejtdQKqTQI5eH6cL6PSUBHKSkIvbUZJAjl4vpwvo9GSj02SREfT7ezhdRCqTQI6etAqcJhudJhPptoiCBHL0pIXsPNnoNHlIIEdBAjl6EsjOk0BOHt2dLiCVSSBHT96ATlMomauXNGQIYhQkkKMnSz8mAwnkZCGBHAUJ5OhVOV2AAJR2ugJhkY1+oyCBHL1KpwsQgFLSRk4O0kKOgnzzoictZIdovBfD3inbI57xvVHea5yuRwCSKVGRb170JJATzHT13x3yzjqjXb0no9Qcp+sR9UiXRRQkkKMngZwAdmt4m7SGk578f4iCBHL05A0YR6ar/56Qd9ZpuzU81+l6RKsuO11AKpNAjp4EcoxpvJfs1nAeyjvW6XpEu0ggR0ECOXoyyiJGTFe/PSHv9aelbzilXXK6gFQmgRw9aSFHwWoNT94W8UzoLa3htHDO6QJSmQRy9CSQO0Baw2nrvNMFpDIJ5OidcbqAVCGt4bSnkRZyVCSQo3fA6QKSndUannVau/pMktZwWisOGEbI6SJSmQRy9AqdLiAZ1baGx+ehMqQ13DkcdrqAVCeBHD1pIddhuvruDXmvP6VdfaRvuPM54nQBqU4COXqFWH1nnXZxG6s1PGlbxDOhFyrjGmCM0zUJR0ggR0lWe4vS40+WlwMnnK7DCaar797KzLs/rvR9RUe8U+fYYSw6r8NOF5DqpIUcG4XAQKeLSARpDYsWSAs5ShLIsXEAmO90EfFk9w2XykgJ0YL9TheQ6iSQYyMtL+xpPJfD3slbpTUs2uACcMjpIlKdBHJspFUgm66++0Le60tkpIRoh+0Bw5B9tKIkgRwbKT8W2W4Nb4t4JvS0W8Ojna5JpJStTheQDiSQY6MQa/fplNstwWoNzzqlXX2vRanZTtcjUtY2pwtIBzLsLQYef7K8AtjkdB1tpfFcDnmnrarwfWV3VZe7R2t3vzkoleN0XSKlSSDHgLSQY+d94Dqni2iJ6epj9w33lZESIpbKgT1OF5EOJJBj533gX5wuoiGrb3jStohnYi4qYxzSNyxib03AMMJOF5EOJJBjZz3WbgnZThcCjVrD0jfcDqFwmHXbt/PpgQOcOXeOiGmSnZXFgL59uX7yZIZcdVW9409KLYgsAAAT7UlEQVSdPcuqzZv57PhxLly+jEspcnNyGDV0KHOmTSM7K6vdNZScOcPHmzdz6NgxLpeX48vMZOiAAcyZNo2Bffs2+ZxjxcW8+/HHnCgpITMjg3EjR7JozhwyvN5Gx5qmyS9feIGKykq++eCDeD1RRcGH0TxZ1JJAjpHHnywPFeT7VgKLnapB47kS9k7aEvFM7Cmt4Y45W1bGs3/+M2fOn6drVhZDBw7E43Zz7sIF9hw8SL+8vHqBvP/wYZ57803CkQi9evRgzLBhhCMRjhUVsXrrVrbv3ctX77uPvNzcNtew5+BBXnjnHcKRCL1zcxnUrx9ny8r49MABdhcWcu9tt3HtmPpDwssuXeLpV14hEokwYsgQyi5eZOPOnZy/cIGH77qr0Wus276dk6dO8cjdd0cbxiCBHDMSyLH1AQ4Esunqsz/kvb5Y+oajUxUK8cyrr3K2rIwFM2eyYMYM3O7agTNXysu5UlFR87mpNa998AHhSISbrruOBTNnopSqOdfzb73FgSNHeHfVKh5csqRNNVy8fJmXli4lHImwcM4c5kydWvPYzn37eGnpUl59/30G9+9PbvfuNY+t2ryZqlCI+xYuZNKYMZimybN//jP7Dx/meHExA/v1qzn2/MWLfLBuHdeOGcPIIUM6/P2qPh2wJdqTCIuMsoit9xP1QhrPFXukxK6qLp8fpd395spIiegs37CBs2VlTB47lptnzaoXxgBZPl+9lu6Zc+cou3QJr8fD/BkzasIYIMPrZcHMmQAcKypqcw1bd+2iKhSyuifqhDHAxNGjuXbMGMKRCKu21M/Ak6dO4XG7mTja+qPI5XIxddw4AI42eP23li/H7XZz+9y5ba6rBSsChmHG4kRCAjmmHn+yfBfQ9v99HWC6+uyvzLzz40rfX4Ui3mlz7K4JEaVwJMLmTz8FYO706W16jsfdtmHnWT5fm+s4XlICwNWDBzf5+IhBgwDYVVh/LtKVigq6ZGbiqvNLwdelC2B9bdV2FRay59AhFs6Z06G+7SZId0UMSZdF7H0APBjLE9p9w1vrjJQYFcvzCzhZUsKVigq6d+tGn549OXLyJPs++4wr5eVkd+3KqCFDGNzgYl6PnBx65+ZSeu4cKzZtYkGdVnJVKMTyDRsAmD5+fJvrqApZOyB1tcO0oepwv3TlChcuXSIn27qGnJuTw+lz57hSUUGW/dzSs2drHgOorKrirRUrGDpgAFOvidlKqW/H6kRCAjkeYhbIVt/wrGLt6iez6OKs+Iy1V22vHj3407JlbNtTf1jt8g0bGDdiBPctXFhzEUwpxb0LF/L711/nw3Xr2L5nD/3y8mou6plac/P113P95MltrqOr3Wo9e+FCk4+fLSur+fe5CxdqAnnMsGHsP3yYt1esYPG8eZy/eJE1W7eSmZHB8IHWyrDvr13L5fJyvvL5z9frXonCpoBhfBaLEwmLBHLsRdWPrPFciXiu3Rr2TsxFZUprOEHK7Yt1h0+cQJsms6dOZcaECWR16cLhEyd4Y/lydhUWkvnRR9xz6601zxvYty/599/PC++8w4mSEs6cP1/z2NWDBjGkf/92hd/wQYPYsXcv2/fs4abrrqs3ZM3UuqZbBaCisrLm39MmTGDbnj1s37uX7Xv31ty/5MYb6ZqVxfGSEtbv2MGCmTPp3bNnzePhcBiXy4XL1aHey5c68iTRPAnkGHv8yfKignzfGuCG9jzPdPU+EPJeXyStYWdobS1UZpom08aPZ9Gc2sEqY6++mm7Z2fzqj39k2549LJg5k572CIfdhYW8vGwZebm5fOWeexjQpw9V4TAHjx5l6apVPP3qq9x98801F9hac+3o0Xy8aRNnzp/nmVdf5fZ58+jbqxdny8p4f+1aik+fxuVyYZpmvf5it8vFo/fdx7Y9e6xxyF4v40aMYPBVV2GaJq998AG9evRg3rRpABw4fJilq1db51OK4YMGseTGG+nVo0d7vm0SyDEmgRwfz9CGQK7TGu6ByhwPjIx/aaIpmRkZNf+e1kSf78C+fbmqb19OlJRw6PhxenbvztmyMl589126ZGbylc9/vuYiWpfMTCaPHUvP7t3535de4p2VK7nm6qtrHm+J1+Ph4bvu4g9vvMHRoiJ+9cILNY+5lOK22bNZuXEjFVVVjc7ncbuZPn58oz7rtdu2UVxayl/dey8ej4fjxcX87o036NOzJw/cfjvllZW8t2YNT73yCt968MF634sWbAgYxtG2HCjaTgI5Pl4Cfgo0eRlbWsPJp0dO7YjBnjlNjx7MzcnhREkJly5fBmDHvn2EIxFGDxvWZNgOueoqeuTkcO7CBU6UlDCijWN+e/Xowd/+5V+y79AhjhQVUVFZSfdu3Rg/ciRZmZksW70al1L06dWr1XOdv3CBD9evZ8q4cQyz+5JXb92KaZp86Y47alrESile++ADduzdy4yJE9tS5ott+mJEu0ggx8HjT5ZfLMj3vUKdi3sad3nEM2mLtIaT01V9+tT8+0pFRc3FtbqulJcD1PTrltkX3rq00KLskplZc872cLtcXDNiBNeMGFHv/m27dwMwqH//JqdEN/TmihV4vV4W1umCKS4tJcvnq9c9Mbh/fwCKTp9uS3kRJJDjQsYhx8+zAKYr70Bl5p0rK32PVoUzps+2w1gkme7Z2QyyZ7MdPHas0ePlFRWcPHUKgAH2WhLd7BEOx4qLmzxneUVF7dCzOrPqOipimqzeaq0DP2vSpFaP//TAAfYeOsTiuXNrhsIBeL1ewuH6awFVD7dr4+XHNwKGcbKNZYt2kECOn+UVXR58r6rLvSO1u/88lIr+f6SIq3kzZgDw0YYNNeEL1mJDr3/0ERVVVVzVp09Na3LciBEorJlwKzduxNS1OxhVVFbyyvvv16xxMaBOCxzg5WXLMH77W9Zt396ojpOnThGpM5kDrHB/6d13KT59mlFDhzJhVMuDbyoqK3lrxQpGDB7caN2Lfnl5VIVCfLK/dk/SrXbLu3+DOpvxi7YcJNpPuizi5PEny3XQ7/8QuLXVg0VSGDt8OLOnTmX1li386sUXGdSvH74uXTheXMzFy5fJyc7m/kWLaoax9cvL46ZZs/hg3TreW7uWzbt20S8vj1A4zPGSEsorKsjMyOC+225rNKys7MIFayKH3Q1S1zsrV1J8+jT9e/cmu2tXrpSXc7SoiKpQiOEDB/LA7be3+rW8t3YtFZWV3HnTTY0emzd9OjvsdTG27dlDeUUFR4uK6Nm9O9eObnU9qn3I7Ly4kUCOr6eA7wOtX14XSWHRnDkM7t+f9du3c7K0lFA4TI9u3bhhyhTmTZvWqG95wcyZDLnqKtbv2MGx4mL2ffYZyl5+c9KYMcyeMqXeBcO2mDR2LDv27ePUmTMcOXmSzIwMBvbrx5SxY5k0dmyr45qPFRezcedObrn++prheXXl5eby5bvv5v21azl49Chut5vxI0eyaO7ctvRL/0o2M40fpbV8b+Mp6Pc/AzzidB1CxMAVYEDAMM63eqToEOlDjr8fA7IalkgHf5Awji8J5DgLGMZe4FWn6xAiSmHgCaeLSHcSyInxI6cLECJKv5eFhOJPAjkBAoaxHXjH6TqE6KAI0qhICAnkxPmh0wUI0UHPBQzjoNNFdAYSyAkSMIx1WGslC5FKpHWcQBLIifUPyIgLkVr+EDCM/a0fJmJBAjmBAoaxA3ja6TqEaKNLwD86XURnIoGceP8MXHS6CCHa4EcBw4jrpr2iPgnkBAsYRgkQdLoOIVpxEDCcLqKzkUB2hgEcdroIIVrweMAwKls/TMSSBLID7Df6t5yuQ4hmvBcwjDfacqBSyqWU+pJS6nWl1AmlVKVS6qxSarNS6gdKqTat59lRSqnDSimtlBoaz9dJFAlkh9hv+D86XYcQDZQD32jLgUqpgcBG4DngDuAo1jIBa4FhWNdLDiqlvtCRQpRSj9hh+2xHnp+KJJCd9Q3gVKtHCZE4323LMDelVE9gFTAVWAGM0FrP0lp/UWt9B9APa4RGFvCCUuqeONacNiSQHRQwjDPA15yuQwjbBtp+Ie/nwFBgE7BIa11vnQutdUhr/QTwONbOUE8rpfJiWGtakkB2WMAwXgFedroO0eldAR4KGEaktQOVUlcD1d0QX9dat7SD60+BT4Ac4G/rnONZuzviEaXUBKXUy0qpYqVURCn1d0qpw8Az9uEP28fqlrowlFK3KKU+VEqVKaWuKKXWK6WWtPB15CmlfqyU2quUKldKXbCf8zWllCObd0ggJ4evA6VOFyE6tW+3Y0beHVjZsUtrvamlA7W1A8bv7E+bCscbsFrZU7C6PpZi/XL4E7DGPuYg8Ns6H6ubOM9fAcuAbKyFvPYCM4HXlFL3NjxYKTUC2Ap8G+gOvAl8DEzAav2/q5TKbOlriwcJ5CQQMIxS4EFkWrVwxlLat3HpVPt2YxuPrw7ta5toeT4K/BcwUmv9gNZ6sdb611rrfwB+Yx+zWmv9SJ2P39DYt4HbtdYz7X7sKcC/YHWX/HsTxz8PDML663SY1voLdt/3GGA/cDPwvTZ+fTEjgZwkAoaxDFnERSTeUeDBdu6T19u+LWnj8dXHuYCeDR7bC/yr1jraxsjPtNZLG9z3H0AZMEIpNbj6TqXUHGA61ozZx+p2uWitjwF/Z3/6daVUQvfDlEBOLt9DdvQViVMJ3BswjNNxfp2WdmV9XWvdar91G7zV8A6tdRVwyP70qjoPzbNv39Ran23iee8CRUA3av8aSAgJ5CQSMAwT+BJw0ulaRKfwzYBhtNgH3IzqAO/bxuOrJ4eYQMMAPNKB12/K0Wbuv2Df1m3pDrBvW9oBpTrIB7RwTMxJICeZgGGcAh7A2sNMiHh5OmAYv+7gc7fYt9e18fgZ9u0OrXXD93V5B2toqD1dHtUt9pa6aVpq1ceNBHISChjGKto4W0qIDtiKNbKno97CCsCxSqnpLR2olFLAQ/anb0bxmrF03L4d3sIxw+zbE3GupR4J5CQVMIxfYV2UECKWjgJLAobR0tjhFmmtC7GGpQH8vJULX98ExmNdQPt5O1+qyr6N9Zjglfbt55RSuQ0fVErdBvTHWg96S8PH40kCObl9B3jB6SJE2jgLLAwYRixafV/HCvfpwDsNF/dRSnmVUv8XKMDqGnhUa93eZQKq6xwbXan1aa1XYQ3F64b1C6VmvLFSagDw3/an/9PKpJeYk0BOYvZQpEewBqwLEY0KrJbxnlicTGt9GpgDbAcWAIVKqTVKqeeVUm8CxcATWH3Ef6G1fqkDL7PePs8Ue/W43yqlfqOU+nIMvoQvYXVdfBE4pJR60a57H9ZY5A+RcciiIXupzruAmPxHEp1SBHggYBhrWj2yHbTWR4FpWJOa3sXqd70XmI213vePsBYd6tCqhlrrSmAh8LZ97r/EmpE3r6XntfHchcBk4D+xuibuBOYDu7CmeC+yXz+hlDWzUSS7oN8/EGtq6dUOlyJSz2MBw3jS6SJE66SFnCIChnEc6zd4ocOliNTyLQnj1CEt5BRjt5SXAyOcrkUkNQ38jYRxapFATkFBv38A9qLgDpcikpMJPBowjGdaPVIkFemySEH2sKX5WKtSCVFXBHhYwjg1SSCnKDuUr8fav0wIsBYL+mLAMP7gdCGiYySQU5i9BdRNyI4jwtrg4MaAYch7IYVJIKc4ewrs/VjjKUXntAeYGTAM+WspxclFvTQS9Pv/BvgZ4Ha6FpEwH2CtaVzmdCEietJCTiMBw/glsIja9WpFevs1sEjCOH1ICzkNBf3+wVj9yjNaO1akpEvA1wKG8XunCxGxJS3kNBQwjKNYC7/81OlaRMztAKZJGKcnaSGnuaDffxfwNNBo3VeRcn4JPB7NWsYiuUkgdwL2dOv/xVo5S6Sec8BfBwzjT60eKVKaBHInEvT7HwYMpLWcSl4GvhEwjBKnCxHxJ4HcyQT9/n5Yf/re5XQtokUnsC7cveF0ISJxJJA7qaDffz/WRb8+rR0rEkoDvwK+EzCMC60dLNKLBHInFvT7uwEB4O+AljaqFImxEeuiXUx39hCpQwJZEPT7h2Dtf/aA07V0UoewfjG+ZO+jKDopCWRRI+j3z8LaJfg6p2vpJM4CPwB+ETCMqtYOFulPAlk0EvT778Bqsc1yupY0dRH4OfDjgGGcd7oYkTwkkEWzgn7/fOCfgJsdLiVdnAF+AvxMglg0RQJZtCro90/HajEvQabbd0Qh1vjvZwOGccXpYkTykkAWbWZf/HsU+ApwlcPlJLsIsBRrhuSbAcMwHa5HpAAJZNFuQb/fDSwG/hprOrasv1zrANbaIb8NGEaR08WI1CKBLKJir5PxF8DngemAcrYiR5wDXgOeCRjGKqeLEalLAlnEjB3Od2OF8xzSu+V8GHjd/lgVMIyws+WIdCCBLOIi6PfnAXcAC4B5wBBnK4paJbAJWAa8HjCMTxyuR6QhCWSREEG/fyhWMM+3b4c5WU8bnAfWAquA1cCmgGFUOluSSHcSyMIRQb+/JzDB/pho344HshNcShhrWNpuYJd9+wmwW6Yxi0STQBZJI+j3K2AwMNS+HYw1vK6//dETyKrz0dyCSGGgAqub4TJQbH8U1fkoBg4C+wOGEYrLFyREO0kgi5RlB3gW4MMK4UqgUsb8ilQlgSyEEElCpsEKIUSSkEAWQogkIYEshBBJQgJZCCGShASyEEIkCQlkIYRIEhLIQgiRJCSQhRAiSUggCyFEkpBAFkKIJCGBLIQQSUICWQghkoQEshBCJAkJZCGESBISyEIIkSQkkIUQIklIIAshRJKQQBZCiCQhgSyEEElCAlkIIZLE/wfhpD59Z9fLywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print (anion_stats)\n",
    "colours=['#991D1D','#8D6608','#857070']\n",
    "matplotlib.rcParams.update({'font.size': 22})\n",
    "plt.pie(anion_stats[0],labels=['Hex','Cubic','Ortho']\n",
    "        ,startangle=90,autopct='%1.1f%%',colors=colours)\n",
    "plt.axis('equal')\n",
    "plt.savefig('Form-perovskites.png')\n",
    "\n",
    "print ('Number of possible charge neutral perovskites from', search[0], 'to', search[len(search)-1], '=', len(charge_balanced))\n",
    "print ('Number of Pauling sensible perovskites from', search[0], 'to', search[len(search)-1], '=', len(pauling_perov))\n",
    "print ('Number of possible cubic perovskites from', search[0], 'to', search[len(search)-1], '=', len(goldschmidt_cubic))\n",
    "print ('Number of possible ortho perovskites from', search[0], 'to', search[len(search)-1], '=', len(goldschmidt_ortho))\n",
    "print ('Number of possible hexagonal perovskites from', search[0], 'to', search[len(search)-1], '=', len(a_too_large))\n",
    "print ('Number of possible non-perovskites from', search[0], 'to', search[len(search)-1], '=', len(A_B_similar))\n",
    "\n",
    "\n",
    "#print goldschmidt_cubic\n",
    "print( \"----------------------------------------------------------------\")\n",
    "print( \"Structures identified with cubic tolerance factor 0.9 < t < 1.0 \")\n",
    "print( \"----------------------------------------------------------------\")\n",
    "for structure in goldschmidt_cubic:\n",
    "    print( structure[0],structure[1],'(HCOO)3')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Approach 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get list of Element objects\n",
    "search = [el for el in smact.ordered_elements(3,87) if \n",
    "                Element(el).oxidation_states]\n",
    "\n",
    "# Covert to list of Species objects\n",
    "all_species = []\n",
    "for el in search:\n",
    "    for oxi_state in Element(el).oxidation_states:\n",
    "        all_species.append(Species(el,oxi_state,\"6_n\"))\n",
    "        \n",
    "# Define lists of interest\n",
    "A_list = [sp for sp in all_species if \n",
    "          (sp.oxidation == -1) and (sp.ionic_radius)]\n",
    "B_list = [sp for sp in all_species if \n",
    "          (4 <= sp.oxidation <= 5) and (sp.ionic_radius)]\n",
    "C_list = [Species('F',-1,4.47)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We define the different categories of list we will populate\n",
    "charge_balanced = []\n",
    "goldschmidt_cubic = []\n",
    "goldschmidt_ortho = []\n",
    "a_too_large = []\n",
    "A_B_similar = []\n",
    "pauling_perov = []\n",
    "anion_stats = []\n",
    "\n",
    "for combo in product(A_list,B_list,C_list):\n",
    "    A, B, C = combo[0], combo[1], combo[2]\n",
    "    # Check for charge neutrality in 1:1:3 ratio\n",
    "    if (1,1,3) in screening.neutral_ratios(\n",
    "        [A.oxidation, B.oxidation, C.oxidation])[1]:\n",
    "        charge_balanced.append(combo)\n",
    "        # Check for pauling test\n",
    "        if screening.pauling_test([A.oxidation, B.oxidation, C.oxidation],\n",
    "                       [A.pauling_eneg, B.pauling_eneg, C.pauling_eneg]):\n",
    "            pauling_perov.append(combo)\n",
    "            # Calculate tolerance factor\n",
    "            tol = (float(A.ionic_radius) + 4.47)/(np.sqrt(2)*(float(B.ionic_radius)+4.47))\n",
    "            if tol > 1.0:\n",
    "                a_too_large.append(combo)\n",
    "            if tol > 0.9 and tol <= 1.0:\n",
    "                goldschmidt_cubic.append([combo])\n",
    "            if tol >= 0.71 and tol < 0.9:\n",
    "                goldschmidt_ortho.append(combo)\n",
    "            if tol < 0.71:\n",
    "                A_B_similar.append(combo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of possible charge neutral perovskites from Li to Fr = 132\n",
      "Number of Pauling senseibe perovskites from Li to Fr = 132\n",
      "Number of possible cubic perovskites from Li to Fr = 40\n",
      "Number of possible ortho perovskites from Li to Fr = 91\n",
      "Number of possible hexagonal perovskites from Li to Fr = 1\n",
      "Number of possible non-perovskites from Li to Fr = 0\n",
      "----------------------------------------------------------------\n",
      "Structures identified with cubic tolerance factor 0.9 < t < 1.0 \n",
      "----------------------------------------------------------------\n",
      "Cl C (HCOO)3\n",
      "Cl Si (HCOO)3\n",
      "Cl S (HCOO)3\n",
      "Br C (HCOO)3\n",
      "Br Si (HCOO)3\n",
      "Br S (HCOO)3\n",
      "Br V (HCOO)3\n",
      "Br Cr (HCOO)3\n",
      "Br Mn (HCOO)3\n",
      "Br Co (HCOO)3\n",
      "Br Ni (HCOO)3\n",
      "Br Ge (HCOO)3\n",
      "Br Se (HCOO)3\n",
      "I Si (HCOO)3\n",
      "I S (HCOO)3\n",
      "I Ti (HCOO)3\n",
      "I V (HCOO)3\n",
      "I Cr (HCOO)3\n",
      "I Mn (HCOO)3\n",
      "I Fe (HCOO)3\n",
      "I Co (HCOO)3\n",
      "I Ni (HCOO)3\n",
      "I Ge (HCOO)3\n",
      "I Se (HCOO)3\n",
      "I Zr (HCOO)3\n",
      "I Nb (HCOO)3\n",
      "I Mo (HCOO)3\n",
      "I Tc (HCOO)3\n",
      "I Ru (HCOO)3\n",
      "I Rh (HCOO)3\n",
      "I Pd (HCOO)3\n",
      "I Sn (HCOO)3\n",
      "I Tb (HCOO)3\n",
      "I Hf (HCOO)3\n",
      "I Ta (HCOO)3\n",
      "I W (HCOO)3\n",
      "I Re (HCOO)3\n",
      "I Os (HCOO)3\n",
      "I Ir (HCOO)3\n",
      "I Pt (HCOO)3\n"
     ]
    }
   ],
   "source": [
    "print ('Number of possible charge neutral perovskites from', search[0], 'to', search[len(search)-1], '=', len(charge_balanced))\n",
    "print ('Number of Pauling sensible perovskites from', search[0], 'to', search[len(search)-1], '=', len(pauling_perov))\n",
    "print ('Number of possible cubic perovskites from', search[0], 'to', search[len(search)-1], '=', len(goldschmidt_cubic))\n",
    "print ('Number of possible ortho perovskites from', search[0], 'to', search[len(search)-1], '=', len(goldschmidt_ortho))\n",
    "print ('Number of possible hexagonal perovskites from', search[0], 'to', search[len(search)-1], '=', len(a_too_large))\n",
    "print ('Number of possible non-perovskites from', search[0], 'to', search[len(search)-1], '=', len(A_B_similar))\n",
    "\n",
    "\n",
    "#print goldschmidt_cubic\n",
    "print( \"----------------------------------------------------------------\")\n",
    "print( \"Structures identified with cubic tolerance factor 0.9 < t < 1.0 \")\n",
    "print( \"----------------------------------------------------------------\")\n",
    "for structure in goldschmidt_cubic:\n",
    "    print( structure[0][0].symbol,structure[0][1].symbol,'(HCOO)3')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
